//emittance
/*
结论：
验证漂移段发射度是不边的，
纵向用共轭坐标z zp计算得到的发射度是不变的，但是phi，W得到的发射度是变化的。
用phi，zp求得的发射度不变，z w空间发射度变化。推测原因是，z，zp是共轭坐标，但是phi w不是。而phi和z是线性关系。
*/
#ifndef EMITTANCE_H
#define EMITTANCE_H
#include <vector>
#include <cmath>
#include <fstream>
#include <sstream>
#include <iostream>
#include <string>
using namespace std;
void emittance(string);
void calEmit(vector<double> *,vector<double> *,double *);
void emittance(string str)
{
    vector<double> x;
    vector<double> xp;
    vector<double> z;
    vector<double> vz;
    vector<double> zp;
    vector<double> emittance;
    double s=1.20090698858;
    double *emit;
    emit=new double;
    ifstream infile(str.c_str());

    stringstream ss;
    string sRow;
    double x1,x2,x3,x4,x5,x6;
    if(!infile)
    {
	cout<<"can't open emittance.cpp input file"<<endl;
	return;
    }
    while(!infile.eof())
    {
	getline(infile,sRow);
	if(sRow=="") continue;
	ss.clear();
	ss.str(sRow);        
	ss>>x1>>x2>>x3>>x4>>x5>>x6;
	x.push_back(x1);
	xp.push_back(x2);
	z.push_back(x5);
	zp.push_back(x6);
    }
    infile.close();
    calEmit(&x,&xp,emit);
    emittance.push_back(*emit);
    calEmit(&z,&zp,emit);
    emittance.push_back(*emit);
    /*
    infile.open("xxp-0.txt");
    if(!infile)
    {
	cout<<"can't open xxp-0 file"<<endl;
	return;
    }
    vector<double> phi;
    vector<double> w;
    vector<double> zp;
    double zAve=0;
    double vzAve=0;
    double gamma;
    while(!infile.eof())
    {
	getline(infile,sRow);
	if(sRow=="") continue;
	ss.clear();
	ss.str(sRow);        
	ss>>x1>>x2>>x3>>x4>>x5>>x6;
	z.push_back(x5);
	vz.push_back(x6);
	zp.push_back(x6/2.9979e8/0.08011655);
	zAve+=x5;
	vzAve+=x6;
        phi.push_back(x5*360/2.9979e8/0.08011655*324e6);
	gamma=1.0/sqrt(1.0-x6*x6/2.9979e8/2.9979e8);
	w.push_back(938.2723*(gamma-1));
    }
    infile.close();
    calEmit(&z,&zp,emit);
    emittance.push_back(*emit);
//    calEmit(&phi,&zp,emit);
//    //calEmit(&phi,&w,emit);
//    emittance.push_back(*emit);
    zAve/=(x.size());
    vzAve/=(x.size());

    for(size_t j=0;j<1;j++)
    {
	srand( (unsigned)time( NULL ) );   
	for(size_t i=0;i<x.size();i++)
	{
	    x[i]=x[i]+xp[i]*s;
            z[i]=z[i]+vz[i]*5.0e-8-zAve-vzAve*5.0e-8;
    //      z[i]=z[i]+zp[i]*s ;
	    phi[i]=z[i]*360/2.9979e8/0.08011655*324e6;
	    //phi[i]=z[i]*360/vz[i]*324e6;
	}
	calEmit(&x,&xp,emit);
	emittance.push_back(*emit);
        calEmit(&z,&zp,emit);
        emittance.push_back(*emit);
//	//calEmit(&phi,&w,emit);
//	calEmit(&phi,&zp,emit);
//	emittance.push_back(*emit);
    }
        */
    cout<<"emittance develop"<<endl;
    cout<<"emittance from "<<str<<endl;
    for(size_t j=0;j<emittance.size();j++)
	cout<<(emittance[j])<<" ";
	//cout<<*(emittance[j])<<" ";
	cout<<endl;
	ofstream of("zzp-0-after.txt");
	for(size_t i=0;i<z.size();i++)
	{
	  of<<z[i]<<" "<<zp[i]<<endl;;
	}
	of.close();
	
}
void calEmit(vector<double> *x,vector<double> *xp,double *emit)
{
    double xAve=0;
    double xpAve=0;
    double xRms=0;
    double xpRms=0;
    double dXXp=0;
    for(size_t i=0;i<x->size();i++)
    {
	xAve+=(*x)[i];
	xpAve+=(*xp)[i];
    }
    xAve/=(x->size());
    xpAve/=(xp->size());
    for(size_t i=0;i<x->size();i++)
    {
	xRms+=((*x)[i]-xAve)*((*x)[i]-xAve);
	xpRms+=((*xp)[i]-xpAve)*((*xp)[i]-xpAve);
	dXXp+=((*x)[i]-xAve)*((*xp)[i]-xpAve);
    }
    xRms/=(x->size());
    xpRms/=(xp->size());
    dXXp/=(xp->size());
    xRms=sqrt(xRms);
    xpRms=sqrt(xpRms);
   *emit=1.0e4*sqrt(xRms*xRms*xpRms*xpRms-dXXp*dXXp);
   double dA,dB;
   dA=-dXXp/(*emit)*1.0e4;
   dB=xRms*xRms/(*emit)*1.0e4;
   cout<<"alph "<<dA<<" dBeta "<<dB<<endl;
}
#endif
